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Abstract 

When analyzing animal movement, it is important to account for interactions be¬ 
tween individuals. However, statistical models for incorporating interaction behavior 
in movement models are limited. We propose an approach that models dependent 
movement by augmenting a dynamic marginal movement model with a spatial point 
process interaction function within a weighted distribution framework. The approach 
is flexible, as marginal movement behavior and interaction behavior can be modeled 
independently. Inference for model parameters is complicated by intractable nor¬ 
malizing constants. We develop a double Metropolis-Hastings algorithm to perform 
Bayesian inference. We illustrate our approach through the analysis of movement 
tracks of guppies {Poecilia reticulata). 

Keywords: auxiliary variable MCMC algorithm, collective motion, biased correlated ran¬ 
dom walk, group navigation, Poecilia reticulata, state-space model 
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1. INTRODUCTION 


Movement models are important for studying animal behavior as they can reveal how ani¬ 
mals use space and interact with the environment. Information on the movement patterns 
of animal species can play an important role in conservation, particularly for migratory 
species (Durban and Pitman 2012). Many methods exist for modeling individual animal 
movement, including models that account for changing behaviors at different locations and 
times by utilizing Markovian switching models (e.g. Harris and Blackwell (2013); Black- 
well (1997)) and models that account for the animal’s preferences for covariates measured 
throughout the territory (e.g. Hooten et al. (2014); Johnson et ah (2013)). 

Interactions between animals can give insight into the structures of animal societies 
(Mersch et al. 2013). Animal species often exhibit herd or school behavior, and even those 
that do not form groups have movement that depends on the behavior of other individ¬ 
uals. Langrock et al. (2014) incorporate dependence by assuming the animals in a herd 
move around a central point, such as a designated group leader or a latent central location. 
Codling and Bode (2014) propose a model that combines individual navigational behavior 
with the tendency to copy the behavior of other nearby individuals by taking a weighted av¬ 
erage of the two behavioral mechanisms. This enables information sharing among neighbors. 
Perna et al. (2014) consider a model that encourages individuals to have a preferential struc¬ 
ture. For example, an individual might tend to stay directly behind another, thus creating 
a leader-follower relationship. Sumpter (2006) gives a broad overview of animal movement, 
including computer simulation models which utilize self propelled particle (SPP) systems 
with specihc movement rules to account for interaction. 
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We propose a model that describes continuous-time dynamics of animal movement 
(Johnson, London, Lea and Durban 2008) while simultaneously allowing for current-location 
based interactions by modeling animal locations as a spatial interacting point process 
(Mpller and Waagepetersen 2004). Point process models allow interaction between ani¬ 
mal locations such as clustering, regularity, or repulsion, through the use of interaction 
functions. This provides a paradigm for modelling different types of interactions between 
animals including collision avoidance, herding behavior, animals that break off into multi¬ 
ple smaller groups, and animals that interact with each other without moving in herds or 
schools. Our model uses a weighted distribution approach to incorporate several features, 
including 

i. directional persistence through a continuous-time biased correlated random walk, 

ii. inter-animal behavior modeled using spatial point process interaction functions, 

iii. observation error in animal locations. 

Other models exist which incorporate one or more of these features; we propose a flexible 
framework for all three. 

[Figure 1 about here.] 

To illustrate our approach we analyze the guppy {Poecilia reticulata) movement data 
of Bode et al. {2012b) in which ten guppies are released in the lower right section of a 
hsh tank, and are attracted to the top left by shelter in the form of shade and rocks. 
A realization of this experiment is shown in Figure [^a) where the interaction between 
guppies is evident, as the guppies remain together in a shoal. To illustrate the need for 
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statistical models incorporating between-animal dependence, Fignre[^b) shows a simulation 
from an independent movement model, as described in Section 2.2. In the simulation, the 
guppies tend to drift apart, so the model does not replicate the shoaling behavior. In Figure 
l^c) we show a simulated realization from our proposed dynamic point process interaction 
(DPPI) model, described in Section 2.4. Each guppy’s marginal movement is modeled as 
a continuous-time biased correlated random walk which results in smooth paths similar 
to the observed guppy paths. Group movement is modeled using the attraction-repulsion 
interaction function of Goldstein et ah (2015). The simulated guppies in Figure [^c) stay 
together in a group, similar to the observed guppies in Figure [^a). 

The rest of this paper is organized as follows. In Section 2, we introduce the general 
modeling framework, and give several examples of point process interaction functions useful 
for modeling group animal movement. In Section 3, we propose a Markov Ghain Monte 
Garlo algorithm to sample from the posterior distributions of model parameters. We de¬ 
scribe a double Metropolis-Hastings algorithm for inference complicated by the intractable 
normalizing function that arises from our point process interaction approach to modeling 
group movement. In Section 4, we examine the performance of our approach by utilizing 
several simulated movement paths. Finally, in Section 5, we use our approach to analyze 
the guppy movement paths of Bode et ah (20126). 

2. MODELING MOVEMENT DYNAMIGS WITH INTERAGTIONS 
In this section, we describe our proposed approach, starting with a continuous-time stochas¬ 
tic model for the dynamics of individual guppy movement. Next, we aggregate the individ¬ 
ual model to incorporate multiple individuals and describe our point process approach to 
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modeling interactions. Finally, we compare our approach to existing methods. 

Let the unobserved states, consisting of the true locations and instantaneous velocities, 
of individuals (l,...,iL) at a given time L be denoted by At^ = ..., and 

let © denote our vector of parameters. We can write an aggregate group movement model 
by assuming independence and multiplying the marginal densities 

k=l 


where f ©) represents a marginal movement model. That is, the individual’s 


(k) (k) 

state at time ti, cx)-. , is modeled conditional on that individual’s state at time fi_i, CKjdj, 
and the k individuals move independently of each other. To model movement interactions, 
we multiply the marginal model by an interaction function, which is a function of the 
pairwise distance between observations at time f*, yielding a joint distribution 


c(©) 


f{Au\Au_,,&) = 


(j) Ak)_ 


where iljjk{cx[^J,cx[’^^-,&) is the interaction function, which we take from methods in point 
process literature. The resulting model is similar to the weighted distribution approach 
to modeling animal movement. Johnson, Thomas, Ver Hoef and Christ (2008) and Lele 
and Keim (2006) utilize this approach to model a resource selection function for animal 
telemetry data which accounts for animals preferentially selecting certain habitats. In our 
method, the animal’s proximity to neighbors, rather than habitat resource covariates, are 
driving movement behavior. Note that c(©) is an intractable normalizing function of 0. 
This complicates posterior evaluation as we will see later. 
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2.1 Marginal Movement Model 


To develop a group movement model with interactions, we start with an existing move¬ 
ment model for an individual, the continuous time biased correlated random walk model 
(CTCRW) from Johnson, London, Lea and Durban (2008). The CTCRW model specihes an 
Ornstein-Uhlenbeck model for velocity, resulting in movement paths that show directional 
persistence, similar to that of the observed guppy movement paths in Figure 1(a). While 
not important for the guppy data, an additional advantage of the CTCRW model is that 
it allows for observations at non-uniform time points. The CTCRW model is flexible, and 
can easily be adjusted to account for complexities in a given data set. For example, Durban 
and Pitman (2012) use the CTCRW model to estimate the displacement velocities of killer 
whalers; Citta et al. (2013) use an adjusted version of the CTCRW model to analyze haul 
out behavior of Eastern Chukchi beluga whales and Kuhn et al. (2014) use the CTCRW 
model to estimate locations of northern fur seals along foraging tracks. 

Let x{t) and y{t) be the observed location coordinates of the animal at time t, 
and /i*'^^(t) be the true unobserved x and y locations of the animal at time t, and 
and the instantaneous x and y directional velocities of the animal at time t. Let s{t) 

be the observed location and cxt the unobserved state at time t, with 

x{t) 

^ yii) y 

We assume that t G and the locations {x(t),y(t)) belong to The x and y elements 
are assumed to be independent, as a positive correlation between x and y velocities, for 
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example, would indicate movement in a northeast or southwest direction. 

To model directional persistence in movement, and are assumed to follow 

independent continuous-time Ornstein-Uhlenbeck processes. We hrst present the CTCRW 
model for one-dimensional movement, focusing on the x coordinate of Equation Our 
development follows that of Johnson, London, Lea and Durban (2008). 

Given a change in time A, the a:-directional velocity is given by 

+ A) = 7i + - 71 ] + 6 (A), (2) 


where '^i(A) is a normal random variable with mean 0 and variance a‘^[l — exp(— 2 /?A)]/ 2 /?, 
represents the variability in the random velocity, 71 describes the directional drift (mean 
velocity) in the x direction, and controls the autocorrelation in velocity. Equation ([^ 
reveals that the updated velocity at time t -|- A -l- A)) is equal to a weighted average 

of the mean drift ( 71 ), and the velocity at time t {v^^\t)) plus a random term with mean 0 . 
Using this parametrization, small values of imply a higher tendency to continue traveling 
with the same velocity over time. The location /i*^^^(t-|-A) is obtained by integrating velocity 
over time 


^t+A 

(t + A) = (t) + {u)du. 

Assuming we have N observations at times {ti, ...,t]s[) , discretization of the continuous 
time model yields the distributions for the unobserved states, 

/ / (.) \ \ 


( (x) \ 

hi: 


,G) 


N 


Ti(/3,A,) 


,^ = 1,...,A, (3) 


+ di(7i,/3,A,),aVi(/3,A,) 

V / V V 

where Aj is the time change between observations i — 1 and i, Ti{/3,Ai) accounts for the 
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directional persistence, 


/ 


Ti(AA,) = 

Aj) models directional drift, 

di(7i,/3, Ai) = 7i 


/3 




0 


A,- 






V 


1 — e 


and the variance matrix of Equation (|^ is given by 

( 

Vi(/9,A.) = 




ni(/?,Ai) n3(/?,Ai) 

yv3{l3,Ai) V2{^,Ai) 


with 


V2(,P,\) 

Mf^^Ai) 


1 - 

2^5 ’ 

1 - 

2^2 • 


Finally, the observed position of the animal is modeled as a Gaussian random variable 
centered at the true location 


7'~ ivM'’,<Ta. 


(x) 2 


where cr|; represents the observation error variance. To aggregate the x and y dimensional 
distributions into a 2-dimensional model, as given in Equation Q, the covariance terms 
between all x and y elements are set to 0. This yields the marginal model for the individual. 













with parameters (/S, 71 , 72 , cr^, ct|;) and distributions 


~ N{ZoLt^,all2) 


( 4 ) 



(5) 



0 0 1 0 y 


For details about the derivation of the model and examples using this model see Johnson, 


London, Lea and Durban (2008). 

2.2 Independent Group Movement Model 

Assuming independent movement between individuals, this model can be easily extended to 
a group setting. For the remainder of the article we assume that the movement parameters 
{(3, 7i, 72, (t|;) are shared by all individuals. 

Assume that we observe K > 1 animals where every individual is observed at each time 
point (fi, ^ 2 ,-••Gv)- The observed locations are denoted by Sti = {s[^\ for 

ti G fi, G, and the unobserved states are denoted At^ = (ckI^, ..., The 

joint distribution for the unobserved states may be expressed as 


N K 



2=1 k=l 


where /3,71,72, cr^) is the density of a normal random variable for the unob¬ 


served state for individual k at time ti, as dehned in Equation ([^. The joint distribution 


for the observed locations conditional on the unobserved states is therefore 


N K 



(7) 


i=l k=l 
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where f a'^) is the density of a normal random variable for the observation error 

for individual k at time ti, as dehned in Equation (|^, 

2.3 Dynamic Point Process Interaction (DPPI) Model 

If we assume independence between individuals, once two animals start to drift apart, there 
is no mechanism to draw the animals back towards each other. To model schooling or herd 
behavior, we propose an approach motivated by spatial point process models. Consider 
Equation ([^, which gives the distribution of the unobserved states of a set of animals at 
the current time point conditional on the locations at the previous time point. To simplify 
notation, let ©i = (/?, 71 , 72 , cx^, (j|;) describe the parameters for the marginal movement 
model, and let ©2 describe the parameters for a spatial point process interaction function 
'0(-). For each pair of locations at the current time point, we multiply the density by a point 
process interaction function 'ijjjk ^5 ; ©2 j which depends only on the pairwise 

Euclidean distance between the current locations, which we dehne to be 6 = 

\/+ {^y '^ — and parameter © 2 . Note that this is not a function of the 

unobserved velocities. Hence we multiply Equation (|^ by the product of our interaction 
functions 

N K 

n n n 

i—lk=2j<k 

which takes values in M+. For two animals i and j, if the value of ipjki') is small, this 
discourages animals from moving to these locations at the same time, similar to a weighted 
distribution approach for resource selection (Johnson, Thomas, Ver Hoef and Christ 2008). 
The ordering of the individuals does not impact the results. 
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The resulting model has joint density given by: 

h ^1) 9 7i, 72, . 

C(©1,©2) ’ 

where h , a^) represents the density of the observed locations conditional on 

the unobserved states from Equation (|^, g |/3, 71 , 72 , a^) represents the density of the 

unobserved states from Equation (|^, represents the interaction function from 

Equation (|^ and c(©i, © 2 ) is the normalizing function required to ensure that the density 
integrates to 1 and is given by the multidimensional integral over the unobserved states: 


c(©i, © 2 ) 



(^e) 9 l/^> 71,72, '4^{At^.,^\&2)dAt^,^ 


The point process interaction function |© 2 ) should be selected based on the assumed 

interaction behavior of the animals being studied. 

Herding or schooling behavior can be generated when individuals repel each other at 
small distances to avoid collisions, attract each other at mid range distances, and behave 
independently when they are a large distance apart. An interaction function that captures 
this behavior is the attraction-repulsion interaction function found in Goldstein et ah (2015). 
This interaction function is given by: 


N K 




2=1 k=2j<k 


with 


'ijj{r;9i,92,03,R) = < 


if 0 < r < i? 


^jJl{r) = 9i - - 92^'^ if R < r < ri 


( 10 ) 


Mr) = 1 + 


(0-i{r-r2))^ 


if 0 > ri 
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Using this parametrization, 6i gives the peak height of the interaction function, 62 gives 


the location of the peak, and 6^ controls the rate at which the function descends after the 


peak. The values ri and r 2 in Equation (10) are the unique real numbers that make '0(r) 


and continuous, given by the solution to the differential equations 


= V’2(n) 

< 


[Figure 2 about here.] 


See Goldstein et ah (2015) for details. Examples of the interaction functions under different 
parameter settings are given in Figure 


2.4 Comparison with Existing Approaches 

There have been several other models proposed to account for interaction behavior in animal 
movement. Let represent the true locations of each of the animals in the group at time 
ti- Gautrais, dost and Theraulaz (2008) utilize a model where the locations at the next 
time point are only dependent on neighbor’s locations at the current time point, so 

the interaction is a function of At.. Since animals generally interact continuously over time 
we prefer a model that allows modeling of group behavior based on the joint distribution 
of the next location of the individuals in the group, resulting in an interaction which is 
a function of At^^^. This results in a reasonable model even if there are long time lags 
between the observations. Additionally, we consider direct estimation of model parameters, 
whereas Gautrais et ah (2008) utilize extensive simulations under different parametrizations 
followed by analysis of group summary statistics. Mann (2011) discuss Bayesian parameter 
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estimation of a SPP system model, where the interaction term in the model is again assumed 
to depend only on the system state at the previous time point. However, the analysis is 
only accurate if the rate of observations matches the rate at which animals update their 
velocity, implicitly assuming that individuals update their velocities at discrete time points 
(McClintock et ah 2014). 

Potts et ah (2014) propose a similar weighted distribution framework that combines 
three different aspects of movement, individual movement, the effect of the environment 
and the interaction with previous behavior of the rest group, to model an individuals next 
location. These factors are modeled by assuming separability and taking the product of 
three different parts, the movement process, the environmental desirability weighting func¬ 
tion, and the collective interaction which can include all information on group movement up 
to and including the recent time point. We extend this framework by using ideas from point 
process statistics to jointly model the probability of members of a group moving to a new 
location rather than considering the group’s recent history. Another recent approach due 
to Langrock et al. (2014) assumes that animals move around a latent centroid to account 
for group dynamics in animal movement. The movement of the individuals is modeled as a 
hidden Markov model with behavioral states. In one state, the animals may be attracted to 
the group centroid, and follow a biased correlated random walk, whereas in an exploratory 
state, the individual might follow a correlated random walk. Instead of the latent centroid 
approach of Langrock et al. (2014), our method deals with the group dynamics by looking 
at the pairwise behavior between the individuals directly, allowing for different types of 
behavior, such as pairs of animals moving together separate from the group. Mann (2011) 
hnds that parameter estimates can be biased if the time lag for the observations does not 
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match the rate at which individuals update their velocities when only the previous locations 
are considered. Our approach does not have this weakness since we model interaction be¬ 
havior dependent on the current joint locations of the group of individuals, rather than just 
the previous locations using point process interaction functions. Johnson et al. (2013) use 
spatio-temporal point process models to study resource selection, but they do not consider 
animal interactions. 

Our weighted distribution approach provides a general approach to modelling movement 
interactions that is not affected by the timescale of the observations due to the joint modeling 
of the locations. This is an improvement over existing methods which model interactions 
based on the most recent locations under a Markovian assumption. In the case of the 
guppies, we are able to model individual movement using existing dynamic models, and 
interaction using existing point process models which provide a natural way of modeling 
the interaction among points in a plane. Both of these types of models have a large literature 
basis and this makes modeling accessible. 

3. MODEL INFERENCE 

Next, we describe a Metropolis-Hastings algorithm to perform Bayesian inference. We select 
priors for each of the parameters that reflect our limited prior information about the model 
parameters. We will use the same priors for both our simulation study and data analysis. 
For 7 i and 72 we specify conjugate normal priors with zero mean and variance equal to 
10"^, 7 r( 7 i) ~ N{0, 10"^) and 7 r( 72 ) ~ N{0, 10^). For the parameters that are restricted to be 
positive we specify truncated normal priors, denoted truncN(/i, with lower bound 
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given hy Bl and density proportional to 


/(a;|/r, Bl) oc exp ) /{x > Bl] 

where I is the indicator function. The priors chosen are given by /3 rNJ truncN(l, 10"^, 0), 

~ truncN(l, 10^, 0) and ~ truncN(l, 10"^, 0). The parameter R was hxed a priori to 
be the minimum distance between individuals across all time points, denoted R. We have 
additional interaction parameters ^i, 62 and 9^. For 9i and 62 we use truncated normal 
priors; 61 ~ truncN(2,10^, 1) and 62 ~ truncN(.R + 1,10"^, R). Finally, since the effect of 9^ 
on the interaction function is minimal for all 9^, greater than one (see Figure we use a 
uniform prior on ( 0 , 1 ) for 9^. 

Inference is straightforward when the point process interactions are not included in 
the model. For the independent group movement model discussed in Section 2.2, we use 
variable-at-a-time Metropolis-Hastings. At each iteration of our MCMC algorithm, we hrst 
update the unobserved states for each individual at each time point, and then each 

of the model parameters (/9, 71 , 72 , cr^, The Kalman hlter can be used for the model 
with no interactions but it can not be easily extended to the general case; thus we focus on 
a more general method for inference. 

We assessed convergence by monitoring Monte Carlo standard errors using the batch 
means procedures, described in Jones et al. (2006) and Flegal et ah (2008), and by comparing 
kernel density estimates of the posterior of the hrst half of the chain and the second half of 
the chain. 

Inference becomes more challenging when interactions are included in the model. With¬ 
out the interaction function ' 0 ( 0 ) the normalizing constant does not depend on the pa- 


15 



rameters, so it can be ignored for Bayesian inference. However, the normalizing function 
in Equation ([^ is a function of all of the model parameters c(©) = c(©i,© 2 ). In the 
Metropolis-Hastings algorithm, using the model likelihood from Equation (|^, and a pro¬ 
posal density we have acceptance probability: 

. A p(©0g(©1©)h(S,,^|A,,^,©0^(A,,^|©0^(A,,^|©0c(©) \ 

Thus, since the normalizing functions do not cancel out we cannot use Metropolis-Hastings 
without accounting for them. 

Many methods have been suggested to deal with this issue in the point process literature, 
however they are often computationally expensive. Besag (1974) proposed an estimation 
method using psuedo-likelihood which does not work well when there is strong interaction. 
Geyer and Thompson (1992) use importance sampling to estimate the normalizing constant, 
however this method only works if the parameter value used in the importance function is 
close to the maximum likelihood estimate of the parameter. Atchade et ah (2013) propose 
an MCMC algorithm for Bayesian inference. Mpller and Waagepetersen (2004) gives an 
overview of several other estimation methods. Here we use the double Metropolis-Hastings 
(MH) algorithm (Liang 2010). This is an approximate version of the auxiliary variable M-H 
algorithm (Mpller et al. 2006; Murray et ah 2012) but avoids perfect sampling (Propp and 
Wilson 1996) which is not possible from our model. The auxiliary variable is approximately 
simulated using a nested MH sampler. This avoids estimation of the normalizing constant 
at the cost of simulating the path using MCMC. The length of the nested MH sampler must 
be large enough so that the distribution of the auxiliary variable is close to that of a perfect 
sampler. 
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The double MH algorithm (Liang 2010) is 


1 . Generate a proposal ©' from some proposal distribution g(©|©') 

2. Generate an auxiliary Y* = {A*^ from a kernel with stationary distribution 

c(©') 

1 ^* is a approximation of a simulated path from the proposal distribution, this is 
accomplished using a MH algorithm. 


3. Accept ©' with probability a = min (1, R{&, ©0); where R{&, ©') is given by 




and iL (©,©', ^ 


, ST ) is the ratio 


H A 


* 

tl:N ’ 



s (KJ^) ^ At:.vlQ) 

A ®') S I®') i> (A;,„ I®') ■ 


In our model, since none of the parameters can be easily separated from the integration 
over the unobserved states; the normalizing function is a function of all model parame¬ 
ters. Thus, we need to use the double Metropolis-Hastings algorithm for each parameter 
update. Therefore, for each parameter update, we use an MH algorithm to simulate a re¬ 
alization of the unobserved states AT and observations ST from our model with the 
proposal parameters, and use this simulation Y* = {AT^ ST^ to estimate the ratio 
H (©, ©', ST^ ^). This requires a simulation of an entire sample path for each new 

proposal parameter. Note that this estimate is only accurate if the value of © is similar to 
the value of ©', so we elect to use variable at a time updates for all parameters, as opposed 
to block updates of ©. 
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Now we consider the DPPI model from Section 2.4 with the attraction-repulsion inter¬ 
action function from Goldstein et al. (2015). In each iteration of our double Metropolis- 
Hastings algorithm we first update the unobserved states, using a four-dimensional 

block Metropolis-Hastings update, where the unobserved state of each hsh j at each time 
point ti, consisting of the true x and y locations and instantaneous velocities, is up¬ 
dated one at a time. Next, we update each each parameter {/3, 71, 72, 6*1, 6*2, 6*3) 

one at a time using a double Metropolis-Hastings update. For each parameter, we use a 
nested MH sampler to generate an auxiliary variable Y* from the DPPI model using the 
current parameters in the MCMC chain and the proposed parameter to be updated. For 
parameters {/3, 71 , 72 , 61 , 62 , 0^) the auxiliary variable is a simulated realization of the 

unobserved states Y* = A*^ ^ and for a‘^ the auxiliary variable also requires a simulated 
realization of the observations Y* = {A^^ ^^). Both of these auxiliary variables are 
generated using a Metropolis-Hastings algorithm. 

The length of the nested MH sampler used to generate the auxiliary variable was de¬ 
termined by examining the distances between the simulated realizations of the observed 
locations ^ as the length is increased. The length was doubled until the average dis¬ 
tance between locations stabilized, resulting in a nested MH sampler of length 200. The 
double Metropolis-Hastings step is time consuming, since it requires a nested Metropolis- 
Hastings sampler for each parameter at each MCMC step. Convergence was determined 
using the same methods as for the independent movement algorithm. 


18 



4. APPLICATION TO SIMULATED DATA 
To test the performance of our double Metropolis-Hastings algorithm we generated sim¬ 
ulated paths from our DPPI model and recovered the true parameters. We simulated 
three group movement paths with starting locations taken from the starting locations 
of the ten guppies in Figure [^a). In all cases the CTCRW movement parameters 0i 
were set to the means of the posterior distributions from Section 5 (/3 = 0.15 ,71 = 
— 1 . 2,72 = 1.5, cr^ = 1.7, (t|; = 0.4). The interaction parameters ©2 were chosen for 
three different scenarios. In scenario 1 (medium interaction), we used the posterior mean 
parameter values from Section 5 = 32,^2^^ = 33 ,^ 3 ^^ = 0.3) to mimic the guppy 

movement. The parameters in scenario 2 were specihed to encourage stronger interaction 
(0^^ = 100 , 6 ^ 2 ^^ = 20 , 6 ^ 3 ^^ = 0.5). The parameters in scenario 3 were specihed to represent 
a weaker interaction ( 6 ^) = 10, 6^2 = 80,6^3 = 0.5). The interaction functions and simu¬ 

lated paths are plotted in Figure The heights of the interaction functions show that the 
second set of parameters (Figure |^b)) results in the strongest interaction, and the third set 
of parameters (Figure [^c)) results in the weakest interaction. In the simulated movement 
paths, it is apparant that Figure |^c) has less interaction, but it is difficult to compare the 
strength of attraction between Figures |^a) and (b) from the plots of the movement paths 
alone. 

[Figure 3 about here.] 

We hrst estimated the parameters using the independent model that assumes that the 
hsh moved independently, as in Section 2.2. The resulting parameter estimates and 95% 
equi-tailed credible intervals are given in Table 
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[Table 1 about here.] 


Our credible intervals for for 71 , 72 , and a'^ include the true parameters for all of the 
simulations. However, in the medium and strong attraction scenarios the credible intervals 
for f3 and cr^ do not contain the truth. This indicates that assuming independence when 
there is actually interaction among the animals can result in biased parameter estimates. 

Next, we used the correct DPPI model to analyze the simulated data. The results are 
given in Table 

[Table 2 about here.] 

From Table we can see that our algorithm accurately recovers the movement parameters 
©1 with the exception of cr^ which falls just outside the 95% credible interval in the strong 
attraction scenario. In Table we are also successful in recovering 6*1 and 62 , but there 
is greater uncertainty in these parameter estimates than in the movement parameters. Al¬ 
though the simulated paths looked similar in Figure we are able to distinguish between 
the medium attraction and the strong attraction scenarios. However the width of the credi¬ 
ble interval increases as attraction increases, indicating it is harder to differentiate between 
levels of attraction as the peak of our attraction-repulsion interaction function increases. 
For 03 , the posterior is very similar to the prior distribution, a uniform distribution on 
( 0 , 1 ), which indicates that there is not enough information in the simulated data to infer 
the parameter. To test the effect that having an incorrect estimate for 9^ would have on 
the other parameter estimates, the double Metropolis-Hastings algorithm was rerun hxing 
03 at several different values (03 = 0 . 05,03 = 0 . 5,03 = 0.9). The resulting posterior distri¬ 
butions for the other parameters remained consistent with our previous results, so the lack 
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of identifiability of 6 ^ does not invalidate our estimates for the other parameters. 

5. GUPPY DATA 

We now use our approach to analyze the guppy shoal data of Bode et ah (20126), available 
online (Bode et ah 2012a), where the individuals show a tendency to interact, as evident 
by the shoaling behavior in Figure [^a). Gravel and shade were added in one corner of 
the tank to attract the guppies, and a group of ten guppies is released in the opposite 
corner. The full trajectories are observed for the guppies from the time they begin moving 
towards the destination until the hrst guppy reaches the target. The guppies were hlmed 
with a standard dehnition camera, recording 10 frames per second, and tracking software 
(SwisTrack; Lochmatter et ah (2008))was used to obtain the coordinates. One realization of 
the experiment is plotted in Figure [^a). The experiment was repeated several times, but we 
focus our analysis on a single realization of the experiment. Bode et ah (20126) calculated a 
summary statistic based on angles of direction to estimate the social interactions of a group. 
A permutation test, which randomly assigned group membership of guppies to artihcial 
experimental trials, found that the social interaction summary statistic was larger in actual 
groups than in artihcially permutated groups in all but 75 out of 10,000 permutations. Bode 
et ah (20126) concluded that the guppies do interact socially. Using our approach, we are 
able to extend the results of Bode et ah (20126) and directly infer parameter values that 
reflect this interaction between hsh. 

We hrst performed inference using the independent movement model from Section 2.2. 
Next we used our double Metropolis-Bastings algorithm to estimate the parameters for the 
DPPI model described in Section 2.3. The priors in both scenarios were selected to be the 
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same as in the simulation example, described in Section 3 . The results are presented in 


Table [H 


[Table 3 about here.] 

The means of the posterior distributions for the the parameters 71 , 72 , and are almost 
identical for the independent and the interaction models. However, the estimates for (3 
and differ slightly. Our results from the simulation study imply that the independent 
model estimates could be inaccurate, since the hsh interact with each other socially (Bode 
et al. 20126). The results for the movement model parameters ©1 indicate that there is 
autocorrelation in the observations over time, the hsh tend to move toward the shelter 
in the upper left corner, and there is appreciable measurement error but it is very small 
in magnitude, since 0.4 pixels is approximately 0.08 cm. This seems reasonable since the 
tracking software used by Bode et al. (20126) is highly accurate (Lochmatter et al. 2008). 

To compare the independent model and the DPPI model, we analyze the distribution of 
pairwise distances from simulated realizations of the two models. In point process statistics, 
Ripley’s K function, which is described in Mpller and Waagepetersen (2004), can be used to 
analyze the attraction or repulsion between points. The K function, however, requires an 
estimate for the intensity of the point process, which does not exist in our model since each 
point has a unique distribution. Instead, we consider the number of pairs of points that lie 
within a distance of d of each other, a monotone function which starts at 0 and ends at the 
total number of pairs of points in the process, dehned by 

N K 

K'(i) = E E E < i) 

i=l k=2 j<k 
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where / represents the indicator function. Larger values of the function indicate that there 
are more pairs of points within that distance of each other; for example larger values of 
K*{d) for small values of d indicate that there is more attraction between points at small 
scales. To test if our htted model is capturing the interaction between guppies, we simulate 
100 movement paths using draws from the posterior densities of the parameters from the 
independent movement model and from the DPPI model. We calculate K*{d) for each of 
the simulated paths, and create 95% pointwise envelopes for the K-functions in the two 
simulation settings by taking the 2.5% and 97.5% quantiles. The K*{d) function is then 
calculated for the data and is compared to the envelopes. The result is plotted in Figure]^ 
The K* (d) function for the guppy data is above the envelope for the independent movement 
model at small distances, indicating that there is more attraction between individuals that 
can be captured in the independent group movement model. When we use the htted DPPI 
model with an attraction-repulsion interaction function, the envelope includes the K*{d) 
function for the guppy data at all distances, indicating that the inclusion of the interaction 
function improves the performance of the model in the case of the guppies. 

[Figure 4 about here.] 

6 . DISCUSSION 

The movement model with point process interactions we have developed allows us to study 
group movement of individuals by considering location-based interactions directly. Our dou¬ 
ble Metropolis-Hastings algorithm for Bayesian inference allows us to accurately estimate 
parameters. We analyze the movement tracks of a shoal of guppies, which was previously 
studied using permutation tests and summary statistics in Bode et al. (20126), and hnd 
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that the DPPI model captures the observed pairwise interactions between guppies. We are 
able to generate paths with similar distributions of pointwise distances between individu¬ 
als using our model, and show that an independent model fails to do so. We have shown 
that ignoring interactions of the guppies from (Bode et ah 20126) leads to unrealistic group 
movement paths and inaccuracies in parameter estimates. 

One drawback of our model is that the simulated paths appear less smooth than the 
actual paths in the data. This could be due to the time-varying behavior of the guppies, 
which is apparent in Figure [^a), as the guppies change direction during their movement. 
Further, the guppies do not all start to move at the same time. Some guppies linger at 
their start location after they have been released. Thus, our assumption of a constant drift 
shared by all hsh may not hold, and including a time-varying drift term that varies across 
individuals in our CTCRW model might better capture the observed movement behavior. 
However, this increased flexibility would exacerbate the computational cost, and without 
incorporating these improvements we are still able to capture the social interactions. 

In future work, we would like to consider the impact of unobserved animals interacting 
with the group. This could potentially result in biased parameter estimates. For example, 
the strength of the attraction to an individual may be overestimated if there are some 
unobserved animals moving in a group, or the range of attraction may be overestimated 
if there are additional unobserved animals between the group members. The locations of 
unobserved animals could be imputed but this would result in additional computational 
difficulties, particularly if the number of unobserved individuals is unknown. 

Analysis on group movement mechanics have focused on three main features: colli¬ 
sion avoidance at small scales, alignment at medium scales, and attraction at larger scales 
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(Gautrais et al. 2008). Our model as presented in Section 2.3 does not explicitly account for 
the alignment behavior. One method to account for the alignment is to model correlation 
between the velocities of different individuals as a function of their pairwise distance at the 
previous time step. Katz et al. (2011), however, hnd that the alignment is automatically in¬ 
duced by the attraction and repulsion behavior, indicating that this might not be necessary 
to add to the model. 

Animal movement models can vary greatly depending on the species being considered. 
In this case, we have only analyzed the movement of guppies, so the results of our analysis 
may not extend directly to other animals with different types of interactions. The flexibility 
to choose a dynamic movement and interaction function provides the potential to model a 
variety of methods of movement, especially when there is prior knowledge of the animal’s 
behavior. 
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Figure 1: Group Movement Paths 


a) Observed Data b) Independent Simulation c) DPPI Simulation 





a) Plotted paths of a shoal of 10 guppies from Bode et ah (20126). 

b) Plotted paths of a simulated realization from the CTCRW model without interactions. 

c) Plotted paths of a simulated realization from the DPPI model with the 
attraction-repulsion point process interaction function. 
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Figure 2: Behavior of Attraction-Repulsion Interaction Function 


a) Height of Peak b) Location of Peak c) Rate of Descent 



Examples of the attraction-repulsion interaction function from Goldstein et ah (2015). 

a) Demonstrates the effect of changing the peak height parameter di. 

b) Demonstrates the effect of changing the peak location parameter 62 . 

c) Demonstrates the effect of changing the rate of descent parameter 9^. 
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Figure 3: Simulated Data under Different Settings 



The attraction-repulsion point process interaction function for the (a)medium, (b)strong, 
and (c)weak simulated realizations of the model; and plots of the simulated paths for the 
(d)medium, (e)strong, and (f)weak interactions. 


34 











Figure 4: Pairwise Distance Envelope 




Estimates of the K*[d) function for the data compared to 95% equi-tailed conhdence 
intervals calculated from simulated paths using parameters drawn from the posterior 
distributions of (a)the CTCRW model assuming no interactions; and (b)The DPPI model 
with the attraction-repulsion interaction function. 
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Table 1: Simulated Model Assuming Independent Movement 


Interaction Strength 

(3 = 0.15 

7i = -1-2 

72 = 1.5 

= 1.7 

al = 0.4 

Medium 

.184 

(.15, .21) 

-1.25 

(-1.57,-0.91) 

1.64 

(1.33,1.98) 

1.94 

(1.76,2.12) 

.389 

(.36, .41) 

Strong 

.210 

(.18, .23) 

-1.11 

(-1.40,-0.83) 

1.30 

(1.02,1.59) 

1.92 

(1.72,2.11) 

.385 

(.35, .41) 

Weak 

.146 

(.12,.16) 

-1.40 

(-1.78,-1.00) 

1.53 

(1.13,1.93) 

1.75 

(1.61,1.93) 

.392 

(.36, .42) 


Posterior means and 95% equi-tailed credible intervals estimated using a variable at a 
time Metropolis-Hastings algorithm assuming there is no interaction between individuals 
on the data simulated from a DPPI model with medium = 32,6*2^^ = 33,6*3^^ = 0.3), 
strong {9^^ = 100,6^2^^ = 20,6*3^^ = 0.5), and weak (6*® = 10,6*2^^ = 80,6'® = 0.5) 
interaction settings. 
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Table 2: Simulated Model Including Interactions 


Interaction Strength 

/9 = 

= 0.15 

71 = - 1-2 

72 = 1-5 

= 1.7 

= 

0.4 

Medinm 


161 

-1.25 

1.64 

1.72 

.404 



(.13,.18) 

(-1.58,-0.90) 

(1.30,2.01) 

(1.57,1.86) 

(.37, .43) 

Strong 


161 

-1.11 

1.32 

1.51 

.413 



(.13,.18) 

(-1.45,-0.79) 

(1.00,1.64) 

(1.38,1.68) 

(.38, .44) 

Weak 


144 

-1.40 

1.53 

1.74 

.391 



(.12,.16) 

(-1.82,-1.01) 

(1.12,1.92) 

(1.59,1.91) 

(.36, .42) 


Interaction Strength 

0 i = 

(32,100,10) 02 

= (33,20,80) 

03 = (0.3, 0.3, 0.5) 



Medinm 



37.5 

33.7 


.408 





(18.1,74.4) 

(29.6,39.1) 

(.050, .954) 



Strong 



66.9 

19.4 


.614 





(31.0,134.2) 

(16.6,21.5) 

(.073, .983) 



Weak 



12.4 

78.7 


.359 





(4.0,33.3) ( 

;20.2,114.4) 

(.011, .947) 



Posterior means and 95% eqni-tailed credible intervals estimated using the donble 
Metropolis-Hastings algorithm on the data simulated from a DPPI model with medium, 
strong, and weak interaction settings. 
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Table 3: Posterior Summary for the Guppy Data 


Model 

/S 

7i 

72 



Indep. 

.159 

-1.18 

1.51 

1.88 

0.384 


(.13,.18) ( 

-1.56,-0.80) 

(1.14,1.89) 

(1.71,2.04) 

(0.35,0.41) 

Interact 

.145 

-1.17 

1.51 

1.75 

0.395 


(.12,.16) ( 

-1.58,-0.77) 

(1.12,1.89) 

(1.60,1.95) 

(0.36,0.42) 


Model 

01 

02 

03 



Interact 

32.0 

32.9 

0.304 




(15.1,58.2) 

(23.4,44.4) 

(0.019,0.921) 



Posterior means and 95% equi-tailed credible intervals for the guppy data of Bode et ah 
{2012b) assuming no interaction and attraction-repulsion point process interactions, 
estimated using variable at a time Metropolis-Bastings and the double 
Metropolis-Bastings algorithm respectively. 
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